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We propose a method to study lattice QCD at finite T and /i. We compare it with direct results and with the 
Glasgow method by using n/=4 QCD at Im(/i)7^0. We locate the critical endpoint (E) of QCD on the Re(/i)-T 
plane. In this study we use nf— 2+1 dynamical staggered quarks with semi-realistic masses on Lt = 4 lattices. 



/. Introduction. — QCD at finite T and fi is of 
fundamental importance, since it describes rele- 
vant features of particle physics in the early uni- 
verse, in neutron stars and in heavy ion collisions. 
Extensive experimental work has been done with 
heavy ion collisions at CERN and Brookhaven 
to explore the ii-T phase diagram. It is a long- 
standing non-perturbative question, whether a 
critical point exists on the fi-T plane, and par- 
ticularly how to tell its location theoretically |^ . 

QCD at finite /i can be formulated on the 
lattice H]; however, standard Monte-Carlo tech- 
niques can not be used. The reason is that for 
Re(/x)^0 the determinant of the Euclidean Dirac 
operator is complex. This fact spoils any impor- 
tance sampling method. 

An attractive approach to alleviate the prob- 
lem is the "Glasgow method" (see e.g. Ref. [|) 
in which the partition function (Z) is expanded 
in powers of exp(/i/T) by using an ensemble of 
configurations weighted by the ^—0 action. Af- 
ter collecting more than 20 million configurations 
only unphysical results were obtained. The rea- 
son is that the fi=0 ensemble does not overlap 
enough with the finite density states of interest. 
A recent review on QCD at finite T-ji is Rcf. 

We propose a method to reduce the overlap 
problem and determine the phase diagram in the 
fi-T plane (for details see |^). We also locate 
the critical point of QCD. (Similar technique was 
successful for determining the endpoint of the hot 
electroweak plasma e.g. on 4D lattices.) 

II. Overlap improving multi-parameter 
reweighting. — Let us study a generic system 
of fermions V' and bosons (f), where the fermion 
Lagrange density is 'ipM{(f>)tp. Integrating over 



the Grassmann fields we get: 

Z{a) = J V(j)cxp[-Sbos{a,(j))] detM(0,a), 

where a denotes a set of parameters of the La- 
grangian. In the case staggered QCD a consists 
of (3, the quark mass (rrig) and ^. For some choice 
of the parameters a=ao importance sampling can 
be done (e.g. for Re(/i)=0). Rewriting the above 
equation one obtains 

Z{a)^ J V(j)exp[~Sbos{ao,(t))] detM(0,ao) 

f r o / .X o / l^^ detM{(j},a) 1 
I detM((/), ao) J 

We treat the curly bracket as an observable - 
which is measured on each configuration- and the 
rest as the measure. It is known that changing 
only one parameter of the ensemble generated at 
ao provides an accurate value for some observ- 
ables only for high statistics. This is ensured 
by rare fluctuations as the mismatched measure 
occasionally sampled the regions where the inte- 
grand is large. This is the so-called overlap prob- 
lem. Having several parameters the set ao can be 
adjusted to ensure a better overlap than obtained 
by varying only one parameter. 

The basic idea of the method as applied to dy- 
namical QCD can be summarized as follows. We 
study the system at Re(/i)=0 around its transi- 
tion point. Using a Glasgow- type technique we 
calculate the determinants for each configuration 
for a set of fi, which, similarly to the Ferrenberg- 
Swendsen method 0, can be used for reweight- 
ing. The average plaquette values can be used 
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Figure 1. The average of ipip at f3— 5.085 as a 
function of Im(/i), for direct results (squares; their 
sizes give the errors), our technique (crosses) and 
Glasgow-type reweighting (dots). 



to perform an additional reweighting in p. Since 
transition configurations were reweighted to tran- 
sition configurations a much better overlap can 
be observed than by reweighting pure hadronic 
configurations to transition ones as done by the 
Glasgow-type techniques (moving along the tran- 
sition line was also suggested by Ref. ]To||). 

///. Illustration and direct test. — We have di- 
rectly tested these ideas in four-flavor QCD with 
mg=0.05 dynamical staggered quarks. 

We first collected 1200 independent Y^AS^ 
configurations at Re(/x)=lm(/i)=0 and some /3 
values and used the Glasgow-reweighting and also 
our technique to study Re(/i)=0, lm{fj,)^0. At 
Re(^)=0, Im(^)7^0 direct simulations are possi- 
ble. After performing these direct simulations as 
well, a clear comparison can be done. Figure 1 
shows the predictions of the three methods for 
the average quark condensates at /3=5.085 as a 
function of Im(^). The predictions of our method 
agree with the direct results, whereas for larger 
Im(/x) the predictions of the Glasgow method are 
by several standard deviations off. Based on these 
experiences we expect that our method can be 
successfully applied at Re{^)y^O. 

IV. rif ^Jf. staggered quarks. — In QCD with 
rif staggered quarks one should change the deter- 
minants to their rif /A power in our two equations. 
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Figure 2. Im(/3^) as a function of the chemical 
potential. 



Importance sampling works also in this case at 
some P and at Re{ii)=0. Since det M is complex, 
one should choose among the Riemann-sheets of 
the fractional power. We solved this problem an- 
alytically ||. 

In the following we keep fi real and look for the 
zeros of Z on the complex (3 plane. At a first or- 
der phase transition the free energy cx log Z{P) is 
non-analytic. Clearly, a phase transition appears 
only in the V— > c» limit, but not in a finite V. 
Nevertheless, Z has zeros at finite V, generating 
the non-analyticity of the free energy, the Lee- 
Yang zeros These are at complex values of 
the parameters, in our case at complex p. For a 
system with a first order transition these zeros ap- 
proach the real axis in the V— > oo limit (detailed 
analysis suggests 1/V scaling). This oo limit 
generates the non-analyticity of the free energy. 
For a system with crossover Z is analytic, and the 
zeros do not approach the real axis in the V— *■ oo 
limit. 

At Tt^O we used Lt=4, Ls=4,6,8 lattices. T=0 
runs were done on 10^- 16 lattices. mu,d—0-0'^5 
and TOs=0.2 were our bare quark masses. 

At T ^ we determined the complex valued 
Lee- Yang zeros, Po, for different V-s as a function 
of /i. Their V— > oo limit was given by a Po{V) = 
P^ + C/V extrapolation. We used 14000, 3600 
and 840 configurations on Ls=4,6 and 8 lattices, 
respectively. Figure 2 shows Im(/3^) as a function 
of 1.1. For small /i-s the extrapolated Im(/3^) is 
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Figure 3. The T-fj, diagram. Direct results 
are given with errorbars. Dotted Une shows the 
crossover, sohd hue the first order transition. The 
box gives the uncertainties of the endpoint. 



We studied the fj,-T phase diagram of QCD 
with dynamical rt/=2+l quarks. Using our 
method we obtained rE=160±3.5 MeV and 
^_B=725±35 MeV for the endpoint. Though fiE 
is too large to be studied at RHIC/LHC, the end- 
point would probably move closer to the /x=0 axis 
when the quark masses get reduced. At /z=0 we 
obtained Tc= 172±3 MeV. Clearly, more work is 
needed to get the final values. One has to extrap- 
olate to zero step-size in the R-algorithm and to 
the thermodynamic, chiral and continuum limits. 
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T22929 /M28413/OM-MU-708 /IKTAll 1 /NIIF. 
This work was in part based on the MILC col- 
laboration's public lattice gauge theory code: 
littp: / /physics. indiana.edu/~sg/milc. htm]. 



inconsistent with a vanishing value, and predicts 
a crossover. Increasing /i the value of Im(/3^) 
decreases, thus the transition becomes consistent 
with a first order phase transition. (Note, that 
systematic overshooting is a finite V effect.) The 
statistical error was determined by jackknife sam- 
ples of the total Ls = 4, 6 and 8 ensembles. Our 
primary result is fiend — 0.375(20). 

To set the physical scale we used a weighted 
average of Rq, rup and a/ct. Note, that (including 
systematics due to finite V) we have (i?o • Wtt) = 
0.73(6), which is at least twice, mu,d is at least 
four times as large as the physical values. 

Figure 3 shows the phase diagram in physical 
units, thus T as a function of /^b, the baryonic 
chemical potential (which is three times larger 
then the quark chemical potential) . The endpoint 
is at Te = 160 ± 3.5 MeV, a^b = 725 ± 35 MeV. 
At HB^O we obtained 172 ± 3 MeV. 

IV. Conclusions. — We proposed a method - 
an overlap improving multi-parameter reweight- 
ing technique- to numerically study non-zero 
and determine the phase diagram in the T-fx 
plane. Our method is applicable to any number 
of Wilson or staggered quarks. As a direct test 
we showed that for lm(/i)^0 the predictions of 
our method are in complete agreement with the 
direct simulations, whereas the Glasgow method 
suffers from the well-known overlap problem. 
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